The nature of the phase transition in dipolar fluids 
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Monte Carlo computer simulations of a quasi two dimensional (2D) dipolar fluid at low and in- 
termediate densities indicate that the structure of the fluid is well described by an ideal mixture of 
self-assembling clusters. A detailed analysis of the topology of the clusters, of their internal energy 
and of their size (or mass) distributions further suggests that the system undergoes a phase transi- 
tion from a dilute phase characterized by a number of disconnected clusters to a condensed phase 
characterized by a network or spanning (macroscopic) cluster that includes most of the particles in 
the system. 

PACS numbers: 61.20.Ja, 61.20.Gy, 75.50.Mm, 68.65.-k 



The condensation of simple fluids results from the free 
energy balance of the high entropy gas and the low en- 
ergy liquid phases. This transition appears to be generic 
in simple fluids interacting through isotropic intermolecu- 
lar potentials that are repulsive at short distances and at- 
tractive otherwise. The dipolar hard sphere fluid (DHS) 
is a model where hard (or soft) spheres with an em- 
bedded central dipole interact through the dipole-dipolc 
potential. As the average dipolar interaction between 
two dipoles (weighted by the Boltzmann factor) is at- 
tractive one may expect a phase behavior analogous to 
that of simple fluids. A recent calculation of the free en- 
ergy of the DHS at several temperatures, using Monte 
Carlo (MC) simulations 0, suggests the presence of an 
isotropic fluid-fluid transition at low densities, lending 
some support to the analogy with simple fluids. How- 
ever, the structure of DHS at those low densities is dras- 
tically different. Numerical simulations of DHS for 
dipolar interaction strengths of the order of the thermal 
energy, have shown that the anisotropy of the dipolar 
potential promotes the formation of self-assembled ag- 
gregates (chains, rings and more complex clusters - see 
fig.n in sharp contrast with the isotropic compact clus- 
ters observed in simple fluids. Moreover, the pair corre- 
lation function of DHS is strongly peaked at contact and 
the internal energy is nearly independent of the density 
at odds with the behavior of simple fluids. 

Association theories 0, 0, IS IE > that include the ef- 
fect of cluster formation in the thermodynamics, describe 
well the slow variation of the internal energy with the 
density and the size (or mass) distribution of the clusters. 
The simplest theoretical approach (based on simulation 
results 0,IE01) assumes that the only effect of the dipo- 
lar interactions is to drive cluster formation. Thus, the 
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FIG. 1: Typical equilibrium configurations of the quasi 2D 



m* = 2.75 and three reduced densities. 



DHS is described as an ideal mixture of self-assembling 
clusters. 

The Hclmholtz free energy density, /, of a system of 
volume V (or area A in 2D), at absolute temperature T 
and with A^ p particles that self assemble to form non- 
interacting clusters of size N (1 < N < N p ), is HQ 



(1) 



JV=1 



where (3 = (/cbT) -1 (fe is Boltzmann's constant). <fi(N) 
and Vq(N) are, respectively, the density and the parti- 
tion function of clusters of size N. The free energy / is 
minimized with respect to the densities <^(iV), subject to 
a normalization condition N p /V = J2n p =i N</>(N), yield- 
ing the set of equations jjj, 



<KjV) = <7(jV)exp (N(3n(N p ,p)) 



(2) 



where fi(N p , p) is the chemical potential of a system with 
N p particles and density p = N p /V. In general [jj], the 



partition function q(N) may be written as 

q(N) = F(N)exp(-0XN), (3) 

where A is the free energy per particle of an infinite clus- 
ter and — \n(VF(N)) is a sublinear correction (in N) 
to that free energy. The substitution of eqs. @ and © 
in the normalization condition results in, 

N p 

p=Y / NF(N)ex P [P( f i(N p ,p)-X)N], (4) 

iV=l 

and defines implicitly p{N p ,p). The phase behavior 
of the system is obtained from the equation of state 
p(p) = limjVp^oo ju(iVp, p), i.e. the thermodynamic limit 
of eq. (@J. This limit depends crucially on the conver- 
gence the sum p t = J2n=i NF(N). If pt diverges, then 
p(p) is an analytic increasing function of p, bounded by A 
(p(p) < A). Consequently, when p t diverges, the system 
does exhibit any phase transition. On the other hand, 
if pt is finite, p(p) converges non-uniformly to the func- 
tion: p(N p = oo,p) (given through eq. 10}) if p < p t ; 
and A if p > pt- Then, for finite pt, fap) is singular at 
p — pt, signaling a phase transition at this density: when 
p < Pt the structure consists of small disconnected clus- 
ters (p(p) < A) and when p > p t any excess particles 
(with respect to pt) condense in an 'infinite', spanning 
cluster (p(p) = A). This transition is similar to a variety 
of other transitions [g - lamellae formation in systems 
of disk-like micelles, emulsification failure in microemul- 
sions, Bosc-Einstein condensation, etc.- where condensa- 
tion is not driven by the interactions between 'particles' 
(or aggregates) . Previous applications of association the- 
ories to the DHS considered chain and ring formation 
only 0, 0, IE an d failed to predict a phase transition. 
Indeed, infinite chains and rings have the same configu- 
rational entropy and internal energy per particle. 
Their behavior is similar to that of self avoiding random 
walks with F(N) oc (TV 7 ^ 1 + N~ 3+a ), where 7 and a 
are universal exponents known from polymer theory 0. 
Since 7 > 1, p% — ► 00 and the absence of a phase transi- 
tion follows from eq. Q). However, chains and rings are 
only the simplest self-assembled clusters in DHS 0, 0] 
and the previous argument may not hold in general. In 
this letter we provide evidence for the occurrence of a 
phase transition in DHSs by analyzing the distribution 
functions of various types of clusters obtained from exten- 
sive Monte Carlo (MC) DHS simulations. Understand- 
ing the existence and nature of this phase transition is 
important for applications based on dispersions of ferro- 
magnetic nanoparticlcs |ll| . where strong dipolar inter- 
actions are present, as well as for theoretical reasons. In 
fact, the interplay between cluster formation and conden- 
sation is a general problem, relevant in a variety of other 
theoretical contexts |l2T |. 

We have performed extensive MC simulations in the 
canonical ensemble for systems of hard spheres with di- 



2 



4 




' ' ' 'kTe 


jVl ' ' ' ' 1 ' ' ' ' ' i = 






Iff 
It)" 4 




~ 1 
/ 

A 




A-^ 200 400 601 

*AA^ 

aAA *\a 


_ 




\ IW . -v - 



100 200 300 400 500 600 



N 

FIG. 2: Mean reduced density of chains (circles, 4>*(N, 2, 0)), 
rings (squares, (j)*(N, 0, 0)) and networks (triangles, <f>^(N)) of 
size N, as obtained from simulations with p* = 0.1. The lines 
are the theoretical predictions for the chains and rings. In- 
set: mean reduced density of networks at p* — 0.05 (inverted 
triangles), p* = 0.0625 (diamonds) and p* =0.1 (triangles). 



amctcr a and dipole strength m, interacting through the 
pair potential, 

o 

TYl 

Udhs = U H s — [3(Ai • ri 2 ){p2 • hi) - Pi ■ fa] ■ (5) 

r 12 

r\i is the distance between spheres 1 and 2, Uhs the 
hard-sphere potential (= oo if ri2 < cr, otherwise), 
f 12 = ''^J' 1 the unit interparticle vector and fa, fa the 
unit vectors in the direction of the dipole moments of 
spheres 1 and 2, respectively. The centers of the spheres 
and their dipole moments were constrained to lie on the 
same plane, and thus the model is a quasi 2D dipolar hard 
sphere fluid . We used N p — 5776 particles and, for the 
isotherm m* = m(a 3 kBT)^ 1 / 2 — 2.75, simulations were 
performed as described in 0, at reduced densities p* = 
cr 2 ^ = 0.025 , 0.03125 , 0.0375 , 0.05 , 0.075, 0.1, 0.15 and 
0.2 (A is the area of the simulation box). Between 100 
(low densities) and 400 (high densities) equilibrium con- 
figurations were generated, resulting from 2 to 5 xlO 7 
MC cycles (i.e. attempts to rotate and move each of the 
N p dipolar spheres). 

Figure^shows snapshots of equilibrium configurations 
and evidence that the structure of this system is non 
trivial. The particles tend to aggregate in clusters of 
several sizes and topologies; some of the clusters exhibit 
linear aggregation only (chains and rings) and all the 
others (networks) have branches even if most of their 
particles are still linearly aggregated. This qualitative 
picture is quantified by defining a circle of diameter 
r c around each particle i: if there are one, two or more 
additional particles within this circle, then i is an end, 
an interior or a junction particle, respectively |l2|. Two 
particles belong to the same cluster if their separation is 
less than r c and the topology of each cluster is determined 
by the number of ends and junctions: rings have interior 
particles only, chains have two ends and no junctions, 



and networks have at least one junction. The cut off r c 
must be ~ a and in this work we took r c = 1.15cr. 

Both the total internal energy of the system and the 
internal energy of the clusters decrease « 0.5% when 
the density is increased from p* = 0.025 to 0.1. At 
all densities, the difference between these two energies 
is of the order of 0.5%, and well defined size distribu- 
tions are observed even though the clusters break and 
recombine during a simulation run. Thus, as in previ- 
ous works |E 0> IE IE j the description of the system as 
an ideal mixture of self-assembling clusters is justified. 
The theory entailed in eqs. ilL'l-SI 11 is easily generalized 
to include TV-clusters with different topologies. Each N- 
cluster is classified according to its number of ends, n±, 
and junctions, n 3 . The partition function of an TV-cluster 
is now Vq(N, n\, 713), with, 

q(N,m,n 3 ) = F(N,n 1 ,n 3 )exp(-Np\(x 1 ,x 3 )), (6) 

where X(x±, x 3 ) is the free energy per particle of an infi- 
nite cluster with a fraction of ends x\ = n\/N and a frac- 
tion of junctions x 3 = n 3 /N, and — /3 _1 lnVF(N, m,n 3 ) 
is the sublinear correction to that free energy. The gen- 
eralization of eq. (J2J) for the density of TV-clusters with 
topology (ni,n 3 ) is, 

<t>(N,m,n 3 ) = F(N,n 1 ,n 3 )exp(pN(ji-\{x 1 ,x 3 ))). (7) 

The partition functions of chains and rings are approx- 
imated by those of self avoiding random walks (SARW) 
with internal energies given by E c h(N) — —Neo + 2ej 
(chains) and E r (N) = -Ne + e r /N (rings) 0. Here, 
eo is a "bond" energy, i.e. the internal energy per par- 
ticle of an infinite chain or ring, e\ is the energy cost 
of an end, and e r /N a long-range correction to the en- 
ergy of a ring. The resulting size distributions are given, 
at fixed m*, by 0(jV, 2, 0) cx TV 7-1 exp(/3(/z - X )N) and 
et>(N,0,Q) cx N- 3+a exp(/3((x - A )iV + /3e r /N), where 
Ao = A(0, 0). In fig. B we plot these predictions (with 
7 = 1.34 and a = 0.5 of 2D SARW0) for p* = 0.1 and 
the simulation results. The remarkable agreement found 
at this density is also observed at p* < 0.1. 

The implementation of this approach requires the 
knowledge of q{N,ni 1 n 3 ) for networks (i.e. for n 3 ^ 0) 
and simulation results for the various cluster distribu- 
tion functions. Unfortunately this is not possible at 
present, as the number of configurations of a network 
of N monomers, with n\ ends and n 3 junctions, is not 
known and simulations do not yield reliable distribution 
functions for clusters with arbitrary topologies. In or- 
der to proceed, we note that an increase in the num- 
ber of junctions and ends of a N cluster increases the 
internal energy, while it decreases the volume available 
to branches, decreasing the translational entropy of the 
network. However, the number of possible branched ar- 
chitectures increases, increasing the configurational en- 
tropy of the network. We may then conjecture that, at 
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FIG. 3: Mean fraction of ends (< m, n > /N, open sym- 
bols) and junctions (< 713 > /N, full symbols) in networks of 
size N, at p* = 0.1 (circles), 0.075 (squares) and 0.0625 (di- 
amonds). Simulation results for clusters with N in the range 
50 — 600. Note that < n\, n > /N may be approximated by 
a power of N, with a negative exponent larger than —1 (cf. 
with the full line 2iV _1 , the fraction of ends in chains of size 
N). 

fixed N, the function q(N,ni,n 3 ) will exhibit a maxi- 
mum for some (ni(N), n 3 (N)). If, in the spirit of the 
saddle-point approximation, only this maximum is con- 
sidered the mean density of networks of size AT, <j> n {N), 
is, 

MN) = F(N, n u n 3 ) exp(/3N(» - A(ii, x 3 ))), (8) 

where x x = fii(N)/N and x 3 = n 3 (N)/N. Substitut- 
ing this (and the corresponding expressions for chains 
and rings) in eq. 10} yields the equation of state. Fur- 
ther analysis requires the knowledge of the functions F, 
fi\ (N) and n 3 (N) . In fig. [3] we plot the mean fraction 
of ends and junctions of iV-networks, at three different 
densities, obtained from the simulations. Approximating 
ni(iV) and n 3 (N) by these quantities, one obtains, in the 
large N limit, n 3 (N) rj C3A, where C3 is constant (if a 
slight dependence on p* is neglected) and fi\ w ciN^, 
with ci and n constant. Comparison with the result for 
chains shows (fig. El that < rj < 1. As a consequence, 
the internal energy of a iV-network can be approximated 
by ~€()N + eiC\N n + e 3 c 3 N, where e 3 is the energy cost 
of a junction. We note that the linear term may be in- 
corporated into the free energy per particle of the infinite 
cluster A(0, C3). The sublinear term must be accounted 
for through the energetic part of F(N). Then, using 
eq. © and the simulation results of fig.|31we may write, 

MN) = F S (N) exp [-/3e lCl N r > + 0N{p - A(0, c 3 ))] , 

(9) 

where (3 1 In VF$ (N) is the excess entropy of a finite ^V- 
network (Fs (N) is expected to scale with some power of 
N 0)- 

Finally we turn to the analysis of the equation of state 
eq. (0} . Using the approximations described previously 
( (</>(N, 2, 0) + <f>(N, 0, 0)) cx (AH-i + N' 3+a ) exp(0(jj, - 
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Ao)) and 6^ given by eq. (O) and replacing the sums by 
integrals eq. J3J is re- written as, 

P= dN(G rc (N) + G n (N))exp[-l3A n {N p ,p)N], 

(10) 

where G rc oc (iV + Af- 2 + a ) exp(-/3A A) and G„ = 
7VF s (7V)exp(-ci/3eiiV). A = A - A(0,c 3 ) is the dif- 
ference between the free energy per particle of an in- 
finite chain (or ring) and that of an infinite network. 
A„ = A(0,c 3 ) — n(N p ,p) is the difference between the 
free energy per particle of an infinite network and that of 
the system. The dependence of these free energies on p* 
is estimated by analysing the simulation results for the 
density distributions of chains and networks (see fig. q). 
According to eq. Q and eq. the slopes of log^ n (JV) 
and lag(f)(N, 2, 0) for large N are proportional to — A„ 
and — A„ — Ao, respectively. These distributions exhibit, 
up to p* = 0.1, the same qualitative behavior seen in fig-n 
the slopes are always negative, indicating that A„ > 
and A„ + Ao > 0; the slope of the chain distribution 
is always more negative than that of networks, meaning 
that A > 0. The dependence of log cj) n ( N) on density 
is depicted in the inset of fig. a the slope increases with 
p* indicating that A„ decreases with increasing density. 
Therefore, one expects A n to vanish at a sufficiency large 
density, while Ao remains positive. Assuming that these 
trens hold for an infinite system, the phase behavior of 
the DHS can be predicted. 

A phase transition occurs when p t = f ± dN(G rc (N) + 
G n (N)) converges. The analysis of the simulation re- 
sults shows that Ao > when A„ = and G n (N) oc 
exp(— ci(3eiN n ) (with < r\ < 1); then p t converges 
and at this density (in the thermodynamic limit, at the 
given m*) the system exhibits a transition to a phase 
where a finite fraction of particles belongs to a macro- 
scopic cluster. The existence of the transition is related 
to the properties of the clusters responsible for the con- 
vergence of pt- Since the internal energy per particle of 
the infinite network (— £0 + ^303) is larger than that of the 
infinite chain (— £o)j Ao > implies that the entropy per 
particle of the infinite network is larger than that of the 
infinite chain. In other words, the difference between the 
entropies must balance the increase in energy due to the 
junctions. On the other hand, the increase in the number 
of ends stabilizes the infinite network, since < 77 < 1. 
Thus, the transition corresponds to the emergence of a 
phase with a large configurational entropy that balances 
the loss of translational entropy and the increase of the 
internal energy. 

Using eq. (|1C)|> we find that, in the thermodynamic 
limit, hm p _ >p -(^) m * is positive, and conclude that the 
transition is discontinuous at m* = 2.75 13]. A second 
order phase transition occurs when the second moment 
of the distribution 4>(N) diverges at pt, as is the case in 
percolation or in Bose-Einstein condensation. 



Recently, the thermodynamics of a self- assembling sys- 
tem of chains and networks (with no rings and no at- 
tractive interactions between clusters) was studied using 
a mean- field (MF) approximation [l2j. The approxima- 
tions used in that work for the distributions of chains and 
networks prevent a quantitative comparison with the re- 
sults of this letter. Nevertheless we stress that the phase 
transition mechanism is the same in both approaches. 

Combining a detailed analysis of extensive MC simu- 
lations and a theory that includes non-linear clusters, we 
confirmed the existence and clarified the nature of the 
phase transition of the low density DHS fluid. We have 
shown that the transition is driven by the formation of 
a macroscopic network, with a large configurational en- 
tropy that overcomes the cost in internal energy and the 
loss of translational entropy of the phase with high con- 
nectivity as predicted by a previous MF calculation 0] • 

The full phase diagram of the DHS at low densities 
may be studied by generalizing the methods described in 
this letter. However, the entropy driven transition may 
be pre-empted by the stabilization of other condensed 
phases (ordered or not) that are not considered in this 
framework. This is unlikely to occur at the low densi- 
ties of the simulations described here but may prevent a 
critical point to be observed. 

We acknowledge D. Levesque and P.I.C. Teixeira for a 
critical reading of this manuscript and for their valuable 
suggestions. 
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